Merge branch 'davidavdav-factfull'

base/linalg/bunchkaufman.jl

@@ -45,3 +45,7 @@ A_ldiv_B!{T<:BlasReal}(B::BunchKaufman{T}, R::StridedVecOrMat{T}) = LAPACK.sytrs
 function A_ldiv_B!{T<:BlasComplex}(B::BunchKaufman{T}, R::StridedVecOrMat{T})
  (issym(B) ? LAPACK.sytrs! : LAPACK.hetrs!)(B.uplo, B.LD, B.ipiv, R)
 end
+
+## reconstruct the original matrix
+## TODO: understand the procedure described at
+## https://www.nag.com/numeric/FL/nagdoc_fl22/pdf/F07/f07mdf.pdf

base/linalg/cholesky.jl

@@ -162,6 +162,7 @@ convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivot
 
 full{T,S}(C::Cholesky{T,S,:U}) = C[:U]'C[:U]
 full{T,S}(C::Cholesky{T,S,:L}) = C[:L]*C[:L]'
+full(F::CholeskyPivoted) = (ip=invperm(F[:p]); (F[:L] * F[:U])[ip,ip])
 
 size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
 size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d)

base/linalg/factorization.jl

@@ -774,3 +774,17 @@ function At_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArr
  TFB = typeof(one(TF)/one(TB))
  At_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
 end
+
+## reconstruct the original matrix
+full(F::QR) = F[:Q] * F[:R]
+full(F::QRCompactWY) = F[:Q] * F[:R]
+full(F::QRPivoted) = (F[:Q] * F[:R])[:,invperm(F[:p])]
+
+## Can we determine the source/result is Real? This is not stored in the type Eigen
+full(F::Eigen) = F.vectors * Diagonal(F.values) / F.vectors
+
+full(F::Hessenberg) = (fq = full(F[:Q]); (fq * F[:H]) * fq')
+
+full(F::Schur) = (F.Z * F.T) * F.Z'
+
+full(F::SVD) = (F.U * Diagonal(F.S)) * F.Vt

base/linalg/ldlt.jl

@@ -57,3 +57,12 @@ function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
  end
  return B
 end
+
+## reconstruct the original matrix, which is tridiagonal
+function full(F::LDLt)
+ e = copy(F.data.ev)
+ d = copy(F.data.dv)
+ e .*= d[1:end-1]
+ d[2:end] += e .* F.data.ev
+ SymTridiagonal(d, e)
+end

base/linalg/lu.jl

@@ -341,3 +341,5 @@ end
 
 /(B::AbstractMatrix,A::LU) = At_ldiv_Bt(A,B).'
 
+## reconstruct the original matrix
+full(F::LU) = (F[:L] * F[:U])[invperm(F[:p]),:]

doc/stdlib/linalg.rst

@@ -39,6 +39,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
  Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, LowerTriangular, UpperTriangular) of A, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example: ``A=factorize(A); x=A\\b; y=A\\C``.
 
+.. function:: full(F)
+
+ Reconstruct the matrix ``A`` from the factorization ``F=factorize(A)``.
+
 .. function:: lu(A) -> L, U, p
 
  Compute the LU factorization of ``A``, such that ``A[p,:] = L*U``.

test/linalg/cholesky.jl

@@ -92,6 +92,7 @@ debug && println("pivoted Choleksy decomposition")
  if isreal(apd)
  @test_approx_eq apd * inv(cpapd) eye(n)
  end
+ @test_approx_eq full(cpapd) apd
  end
  end
 end

test/linalg/symmetric.jl

@@ -36,6 +36,7 @@ let n=10
  eig(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) # same result, but checks that method works
  @test_approx_eq eigvals(Hermitian(asym), 1:2) d[1:2]
  @test_approx_eq eigvals(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2])) d[1:2]
+ @test_approx_eq full(eigfact(asym)) asym
 
  # relation to svdvals
  @test sum(sort(abs(eigvals(Hermitian(asym))))) == sum(sort(svdvals(Hermitian(asym))))
@@ -81,4 +82,3 @@ let A = [1.0+im 2.0; 2.0 0.0]
  @test !ishermitian(A)
  @test_throws ArgumentError Hermitian(A)
 end
-

test/linalg1.jl

@@ -61,6 +61,7 @@ debug && println("(Automatic) Square LU decomposition")
  @test_approx_eq (l*u)[invperm(p),:] a
  @test_approx_eq a * inv(lua) eye(n)
  @test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
+ @test_approx_eq full(lua) a
 
 debug && println("Thin LU")
  lua = lufact(a[:,1:n1])
@@ -77,6 +78,7 @@ debug && println("QR decomposition (without pivoting)")
  @test_approx_eq q*full(q, thin=false)' eye(n)
  @test_approx_eq q*r a
  @test_approx_eq_eps a*(qra\b) b 3000ε
+ @test_approx_eq full(qra) a
 
 debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
  qrpa = factorize(a[1:n1,:])
@@ -87,6 +89,7 @@ debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is onl
  @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:n1,p] : a[1:n1,:]
  @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:n1,:]
  @test_approx_eq_eps a[1:n1,:]*(qrpa\b[1:n1]) b[1:n1] 5000ε
+ @test_approx_eq full(qrpa) a[1:5,:]
 
 debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats
  qrpa = factorize(a[:,1:n1])
@@ -96,6 +99,7 @@ debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is on
  @test_approx_eq q*full(q, thin=false)' eye(n)
  @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:n1]
  @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:n1]
+ @test_approx_eq full(qrpa) a[:,1:5]
 
 debug && println("non-symmetric eigen decomposition")
  if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
@@ -132,6 +136,7 @@ debug && println("Schur")
  @test_approx_eq sort(real(f[:values])) sort(real(d))
  @test_approx_eq sort(imag(f[:values])) sort(imag(d))
  @test istriu(f[:Schur]) || iseltype(a,Real)
+ @test_approx_eq full(f) a
  end
 
 debug && println("Reorder Schur")
@@ -177,6 +182,7 @@ debug && println("singular value decomposition")
  if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
  usv = svdfact(a)
  @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a
+ @test_approx_eq full(usv) a
  end
 
 debug && println("Generalized svd")

test/linalg2.jl

@@ -69,6 +69,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
  Tldlt = ldltfact(Ts)
  x = Tldlt\v
  @test_approx_eq x invFsv
+ @test_approx_eq full(full(Tldlt)) Fs
  end
 
  # eigenvalues/eigenvectors of symmetric tridiagonal