Fix JuliaLang#10970. Allow LU factorizing non-square sparse matrices.

pwl · Apr 27, 2015 · d354afc · d354afc
1 parent 982c60a
commit d354afc
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 9 deletions.
diff --git a/base/sparse/umfpack.jl b/base/sparse/umfpack.jl
@@ -102,7 +102,6 @@ type UmfpackLU{Tv<:UMFVTypes,Ti<:UMFITypes} <: Factorization{Tv}
 end
 
 function lufact{Tv<:UMFVTypes,Ti<:UMFITypes}(S::SparseMatrixCSC{Tv,Ti})
- S.m == S.n || error("argument matrix must be square")
 
  zerobased = S.colptr[1] == 0
  res = UmfpackLU(C_NULL, C_NULL, S.m, S.n,
@@ -270,8 +269,8 @@ for itype in UmfpackIndexTypes
  end
  function umf_extract(lu::UmfpackLU{Float64,$itype})
  umfpack_numeric!(lu) # ensure the numeric decomposition exists
- (lnz,unz,n_row,n_col,nz_diag) = umf_lunz(lu)
- Lp = Array($itype, n_col + 1)
+ (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+ Lp = Array($itype, n_row + 1)
  Lj = Array($itype, lnz) # L is returned in CSR (compressed sparse row) format
  Lx = Array(Float64, lnz)
  Up = Array($itype, n_col + 1)
@@ -289,14 +288,14 @@ for itype in UmfpackIndexTypes
  Up,Ui,Ux,
  P, Q, C_NULL,
  &0, Rs, lu.numeric)
- (transpose(SparseMatrixCSC(n_row,n_row,increment!(Lp),increment!(Lj),Lx)),
- SparseMatrixCSC(n_row,n_col,increment!(Up),increment!(Ui),Ux),
+ (transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), Lx)),
+ SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), Ux),
  increment!(P), increment!(Q), Rs)
  end
  function umf_extract(lu::UmfpackLU{Complex128,$itype})
  umfpack_numeric!(lu) # ensure the numeric decomposition exists
- (lnz,unz,n_row,n_col,nz_diag) = umf_lunz(lu)
- Lp = Array($itype, n_col + 1)
+ (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+ Lp = Array($itype, n_row + 1)
  Lj = Array($itype, lnz) # L is returned in CSR (compressed sparse row) format
  Lx = Array(Float64, lnz)
  Lz = Array(Float64, lnz)
@@ -316,8 +315,8 @@ for itype in UmfpackIndexTypes
  Up,Ui,Ux,Uz,
  P, Q, C_NULL, C_NULL,
  &0, Rs, lu.numeric)
- (transpose(SparseMatrixCSC(n_row,n_row,increment!(Lp),increment!(Lj),complex(Lx,Lz))),
- SparseMatrixCSC(n_row,n_col,increment!(Up),increment!(Ui),complex(Ux,Uz)),
+ (transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), complex(Lx, Lz))),
+ SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), complex(Ux, Uz)),
  increment!(P), increment!(Q), Rs)
  end
  end

diff --git a/test/sparsedir/umfpack.jl b/test/sparsedir/umfpack.jl
@@ -42,6 +42,15 @@ for Ti in Base.SparseMatrix.UMFPACK.UMFITypes.types
  @test_approx_eq scale(Rs,Ac)[p,q] L*U
 end
 
+for elty in (Float64, Complex128)
+ for (m, n) in ((10,5), (5, 10))
+ A = sparse([1:min(m,n); rand(1:m, 10)], [1:min(m,n); rand(1:n, 10)], elty == Float64 ? randn(min(m, n) + 10) : complex(randn(min(m, n) + 10), randn(min(m, n) + 10)))
+ F = lufact(A)
+ L, U, p, q, Rs = F[:(:)]
+ @test_approx_eq scale(Rs,A)[p,q] L*U
+ end
+end
+
 #4523 - complex sparse \
 x = speye(2) + im * speye(2)
 @test_approx_eq ((lufact(x) \ ones(2)) * x) (complex(ones(2)))