diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md
index 52e7860999287..d1e7629566354 100644
--- a/stdlib/LinearAlgebra/docs/src/index.md
+++ b/stdlib/LinearAlgebra/docs/src/index.md
@@ -394,6 +394,7 @@ LinearAlgebra.diagind
 LinearAlgebra.diag
 LinearAlgebra.diagm
 LinearAlgebra.rank
+LinearAlgebra.rref
 LinearAlgebra.norm
 LinearAlgebra.opnorm
 LinearAlgebra.normalize!
diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index 4eee4ecd92b11..8c573b826b921 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -126,6 +126,8 @@ export
     rdiv!,
     reflect!,
     rotate!,
+    rref,
+    rref!,
     schur,
     schur!,
     svd,
@@ -194,7 +196,7 @@ julia> LinearAlgebra.stride1(B)
 2
 ```
 """
-stride1(x) = stride(x,1)
+stride1(x) = stride(x, 1)
 stride1(x::Array) = 1
 stride1(x::DenseArray) = stride(x, 1)::Int
 
@@ -219,7 +221,7 @@ julia> LinearAlgebra.checksquare(A, B)
 ```
 """
 function checksquare(A)
-    m,n = size(A)
+    m, n = size(A)
     m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
     m
 end
@@ -227,8 +229,8 @@ end
 function checksquare(A...)
     sizes = Int[]
     for a in A
-        size(a,1)==size(a,2) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(a))"))
-        push!(sizes, size(a,1))
+        size(a, 1) == size(a, 2) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(a))"))
+        push!(sizes, size(a, 1))
     end
     return sizes
 end
@@ -375,6 +377,7 @@ include("bunchkaufman.jl")
 include("diagonal.jl")
 include("bidiag.jl")
 include("uniformscaling.jl")
+include("rref.jl")
 include("hessenberg.jl")
 include("givens.jl")
 include("special.jl")
@@ -407,10 +410,10 @@ of the problem that is solved on each processor.
     the standard library `InteractiveUtils`.
 """
 function peakflops(n::Integer=2000; parallel::Bool=false)
-    a = fill(1.,100,100)
-    t = @elapsed a2 = a*a
-    a = fill(1.,n,n)
-    t = @elapsed a2 = a*a
+    a = fill(1., 100, 100)
+    t = @elapsed a2 = a * a
+    a = fill(1., n, n)
+    t = @elapsed a2 = a * a
     @assert a2[1,1] == n
     if parallel
         let Distributed = Base.require(Base.PkgId(
@@ -418,8 +421,8 @@ function peakflops(n::Integer=2000; parallel::Bool=false)
             return sum(Distributed.pmap(peakflops, fill(n, Distributed.nworkers())))
         end
     else
-        return 2*Float64(n)^3 / t
-    end
+        return 2 * Float64(n)^3 / t
+end
 end
 
 
@@ -428,9 +431,9 @@ function versioninfo(io::IO=stdout)
         openblas_config = BLAS.openblas_get_config()
         println(io, "BLAS: libopenblas (", openblas_config, ")")
     else
-        println(io, "BLAS: ",Base.libblas_name)
+        println(io, "BLAS: ", Base.libblas_name)
     end
-    println(io, "LAPACK: ",Base.liblapack_name)
+    println(io, "LAPACK: ", Base.liblapack_name)
 end
 
 function __init__()
diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl
index 8e0ca4fb72ad5..0c439ca8332a7 100644
--- a/stdlib/LinearAlgebra/src/generic.jl
+++ b/stdlib/LinearAlgebra/src/generic.jl
@@ -1011,6 +1011,78 @@ function rank(A::AbstractMatrix; atol::Real = 0.0, rtol::Real = (min(size(A)...)
 end
 rank(x::Number) = x == 0 ? 0 : 1
 
+"""
+    rref(A)
+Compute the reduced row echelon form of the matrix A.
+Since this algorithm is sensitive to numerical imprecision,
+* Complex numbers are converted to ComplexF64
+* Integer, Float16 and Float32 numbers are converted to Float64
+* Rational are kept unchanged
+```jldoctest
+julia> rref([ 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
+julia> rref([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
+julia> rref([ 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
+```
+"""
+function rref!(A::Matrix{T}, ɛ=T <: Union{Rational,Integer} ? 0 : eps(norm(A, Inf))) where T
+    nr, nc = size(A)
+    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] .= zero(T)
+            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
+            i += 1
+            j += 1
+        end
+    end
+    return A
+end
+
+rrefconv(::Type{T}, A::Matrix) where {T} = rref!(copyto!(similar(A, T), A))
+rref(A::Matrix{T}) where {T} = rref!(copy(A))
+rref(A::Matrix{T}) where {T <: Complex} = rrefconv(ComplexF64, A)
+rref(A::Matrix{ComplexF64}) = rref!(copy(A))
+rref(A::Matrix{T}) where {T <: Union{Integer,Float16,Float32}} = rrefconv(Float64, A)
+rref(A::AbstractMatrix) = rref(Matrix(A))
+
 """
     tr(M)