diff --git a/base/linalg/bunchkaufman.jl b/base/linalg/bunchkaufman.jl
index 891c54768c7d5..84613c3b5d628 100644
--- a/base/linalg/bunchkaufman.jl
+++ b/base/linalg/bunchkaufman.jl
@@ -43,3 +43,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
+## http://www.nag.com/numeric/FL/nagdoc_fl22/pdf/F07/f07mdf.pdf
diff --git a/base/linalg/cholesky.jl b/base/linalg/cholesky.jl
index 8c71fca854567..4ae3624313c0a 100644
--- a/base/linalg/cholesky.jl
+++ b/base/linalg/cholesky.jl
@@ -108,6 +108,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)
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index e38d9b973d044..cd50a16e97231 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -716,3 +716,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
diff --git a/base/linalg/ldlt.jl b/base/linalg/ldlt.jl
index 09faba588e108..863cf2bd6f602 100644
--- a/base/linalg/ldlt.jl
+++ b/base/linalg/ldlt.jl
@@ -57,3 +57,6 @@ function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
     end
     return B
 end
+
+## reconstruct the original matrix, but in full rather than tridiagonal
+full(F::LDLt) = (L=Triangular(F.data, :L, true); L * Diagonal(diag(F.data)) * L')
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
index f707b697cddbd..b5e31184d70e7 100644
--- a/base/linalg/lu.jl
+++ b/base/linalg/lu.jl
@@ -318,3 +318,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]),:]
diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst
index 5bec267b1f210..3f6106ed7f13f 100644
--- a/doc/stdlib/linalg.rst
+++ b/doc/stdlib/linalg.rst
@@ -41,6 +41,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    ``factorize!`` is the same as :func:`factorize`, but saves space by overwriting the input ``A``, instead of creating a copy.
 
+.. 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``.
diff --git a/test/linalg1.jl b/test/linalg1.jl
index 31057b0d3ca93..820834342971a 100644
--- a/test/linalg1.jl
+++ b/test/linalg1.jl
@@ -78,6 +78,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
 
 debug && println("(Automatic) Bunch-Kaufman factor of indefinite matrix")
@@ -100,6 +101,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:5])
@@ -116,6 +118,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:5,:])
@@ -126,6 +129,7 @@ debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is onl
     @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:5,p] : a[1:5,:]
     @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:5,:]
     @test_approx_eq_eps a[1:5,:]*(qrpa\b[1:5]) b[1:5] 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:5])
@@ -135,6 +139,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:5]
     @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:5]
+    @test_approx_eq full(qrpa) a[:,1:5]
 
 debug && println("symmetric eigen-decomposition")
     if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
@@ -146,6 +151,7 @@ debug && println("symmetric eigen-decomposition")
         @test_approx_eq abs(eigfact(Hermitian(asym), d[1]-10*eps(d[1]), d[2]+10*eps(d[2]))[:vectors]'v[:,1:2]) eye(eltya, 2)
         @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
     end
 
 debug && println("non-symmetric eigen decomposition")
@@ -178,6 +184,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")
@@ -205,6 +212,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")
diff --git a/test/linalg2.jl b/test/linalg2.jl
index f9cff896fa847..995427d2a7aa2 100644
--- a/test/linalg2.jl
+++ b/test/linalg2.jl
@@ -82,6 +82,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
         Tldlt = ldltfact(Ts)
         x = Tldlt\v
         @test_approx_eq x invFsv
+        @test_approx_eq full(Tldlt) Fs
     end
 
     # eigenvalues/eigenvectors of symmetric tridiagonal