From aae89599f3b2e1bcd906d314ef879b0fa6e5d25b Mon Sep 17 00:00:00 2001
From: crstnbr <crstnbr@gmail.com>
Date: Wed, 13 Feb 2019 10:46:07 +0100
Subject: [PATCH 1/6] alg keyword for LinearAlgebra.svd

---
 NEWS.md                                |  1 +
 stdlib/LinearAlgebra/src/bidiag.jl     |  4 +--
 stdlib/LinearAlgebra/src/svd.jl        | 42 ++++++++++++++++++--------
 stdlib/LinearAlgebra/src/triangular.jl |  2 +-
 stdlib/LinearAlgebra/test/svd.jl       | 31 +++++++++++++++++++
 5 files changed, 64 insertions(+), 16 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 2ed6586a110dc..303864055e955 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -45,6 +45,7 @@ Standard library changes
 * The BLAS submodule no longer exports `dot`, which conflicts with that in LinearAlgebra ([#31838]).
 * `diagm` and `spdiagm` now accept optional `m,n` initial arguments to specify a size ([#31654]).
 * `Hessenberg` factorizations `H` now support efficient shifted solves `(H+µI) \ b` and determinants, and use a specialized tridiagonal factorization for Hermitian matrices. There is also a new `UpperHessenberg` matrix type ([#31853]).
+* Added keyword argument `alg` to `svd` and `svd!` that allows one to switch between different SVD algorithms ([#31057]).
 
 #### SparseArrays
 
diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl
index 8547e1c2e5f0f..7f409c95c7063 100644
--- a/stdlib/LinearAlgebra/src/bidiag.jl
+++ b/stdlib/LinearAlgebra/src/bidiag.jl
@@ -198,8 +198,8 @@ function svd!(M::Bidiagonal{<:BlasReal}; full::Bool = false)
     d, e, U, Vt, Q, iQ = LAPACK.bdsdc!(M.uplo, 'I', M.dv, M.ev)
     SVD(U, d, Vt)
 end
-function svd(M::Bidiagonal; full::Bool = false)
-    svd!(copy(M), full = full)
+function svd(M::Bidiagonal; kw...)
+    svd!(copy(M), kw...)
 end
 
 ####################
diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 89b7501811230..7a33e31d57cbc 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -54,6 +54,12 @@ function SVD{T}(U::AbstractArray, S::AbstractVector{Tr}, Vt::AbstractArray) wher
         convert(AbstractArray{T}, Vt))
 end
 
+
+abstract type SVDAlgorithm end
+struct DivideAndConquer <: SVDAlgorithm end
+struct Simple <: SVDAlgorithm end
+
+
 # iteration for destructuring into components
 Base.iterate(S::SVD) = (S.U, Val(:S))
 Base.iterate(S::SVD, ::Val{:S}) = (S.S, Val(:V))
@@ -61,7 +67,7 @@ Base.iterate(S::SVD, ::Val{:V}) = (S.V, Val(:done))
 Base.iterate(S::SVD, ::Val{:done}) = nothing
 
 """
-    svd!(A; full::Bool = false) -> SVD
+    svd!(A; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) -> SVD
 
 `svd!` is the same as [`svd`](@ref), but saves space by
 overwriting the input `A`, instead of creating a copy.
@@ -92,18 +98,25 @@ julia> A
   0.0       0.0  -2.0  0.0  0.0
 ```
 """
-function svd!(A::StridedMatrix{T}; full::Bool = false) where T<:BlasFloat
+function svd!(A::StridedMatrix{T}; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) where T<:BlasFloat
     m,n = size(A)
     if m == 0 || n == 0
         u,s,vt = (Matrix{T}(I, m, full ? m : n), real(zeros(T,0)), Matrix{T}(I, n, n))
     else
-        u,s,vt = LAPACK.gesdd!(full ? 'A' : 'S', A)
+        if typeof(alg) == DivideAndConquer
+            u,s,vt = LAPACK.gesdd!(full ? 'A' : 'S', A)
+        elseif typeof(alg) == Simple
+            c = full ? 'A' : 'S'
+            u,s,vt = LAPACK.gesvd!(c, c, A)
+        else
+            throw(ArgumentError("Unsupported value for `alg` keyword."))
+        end
     end
     SVD(u,s,vt)
 end
 
 """
-    svd(A; full::Bool = false) -> SVD
+    svd(A; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) -> SVD
 
 Compute the singular value decomposition (SVD) of `A` and return an `SVD` object.
 
@@ -120,6 +133,9 @@ and `V` is `N \\times N`, while in the thin factorization `U` is `M
 \\times K` and `V` is `N \\times K`, where `K = \\min(M,N)` is the
 number of singular values.
 
+If `alg = DivideAndConquer()` (default) a divide-and-conquer algorithm is used to calculate the SVD.
+One can set `alg = Simple()` to use a simple (typically slower) algorithm instead.
+
 # Examples
 ```jldoctest
 julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
@@ -144,21 +160,21 @@ julia> u == F.U && s == F.S && v == F.V
 true
 ```
 """
-function svd(A::StridedVecOrMat{T}; full::Bool = false) where T
-    svd!(copy_oftype(A, eigtype(T)), full = full)
+function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) where T
+    svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
 end
-function svd(x::Number; full::Bool = false)
+function svd(x::Number; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
     SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
 end
-function svd(x::Integer; full::Bool = false)
-    svd(float(x), full = full)
+function svd(x::Integer; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+    svd(float(x), full = full, alg = alg)
 end
-function svd(A::Adjoint; full::Bool = false)
-    s = svd(A.parent, full = full)
+function svd(A::Adjoint; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+    s = svd(A.parent, full = full, alg = alg)
     return SVD(s.Vt', s.S, s.U')
 end
-function svd(A::Transpose; full::Bool = false)
-    s = svd(A.parent, full = full)
+function svd(A::Transpose; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+    s = svd(A.parent, full = full, alg = alg)
     return SVD(transpose(s.Vt), s.S, transpose(s.U))
 end
 
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
index fc78bf57547ed..e751a4c801d5c 100644
--- a/stdlib/LinearAlgebra/src/triangular.jl
+++ b/stdlib/LinearAlgebra/src/triangular.jl
@@ -2513,7 +2513,7 @@ eigen(A::AbstractTriangular) = Eigen(eigvals(A), eigvecs(A))
 # Generic singular systems
 for func in (:svd, :svd!, :svdvals)
     @eval begin
-        ($func)(A::AbstractTriangular) = ($func)(copyto!(similar(parent(A)), A))
+        ($func)(A::AbstractTriangular; kwargs...) = ($func)(copyto!(similar(parent(A)), A); kwargs...)
     end
 end
 
diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl
index 7a29d86974a61..0a23720bd6cc5 100644
--- a/stdlib/LinearAlgebra/test/svd.jl
+++ b/stdlib/LinearAlgebra/test/svd.jl
@@ -139,4 +139,35 @@ aimg  = randn(n,n)/2
     end
 end
 
+
+
+@testset "SVD Algorithms" begin
+    ≊(x,y) = isapprox(x,y,rtol=1e-15)
+
+    allpos = (v) -> begin
+        for e in v
+            e < 0 && return false
+        end
+        return true
+    end
+
+    x = [0.1 0.2; 0.3 0.4]
+
+    for alg in [LinearAlgebra.Simple(), LinearAlgebra.DivideAndConquer()]
+        sx1 = svd(x, alg = alg)
+        @test sx1.U * Diagonal(sx1.S) * sx1.Vt ≊ x
+        @test sx1.V * sx1.Vt ≊ I
+        @test sx1.U * sx1.U' ≊ I
+        @test allpos(sx1.S)
+
+        sx2 = svd!(copy(x), alg = alg)
+        @test sx2.U * Diagonal(sx2.S) * sx2.Vt ≊ x
+        @test sx2.V * sx2.Vt ≊ I
+        @test sx2.U * sx2.U' ≊ I
+        @test allpos(sx2.S)
+    end
+end
+
+
+
 end # module TestSVD

From bbbff75d4a70233fcfc76be5f38cc94f967c513b Mon Sep 17 00:00:00 2001
From: Carsten Bauer <crstnbr@gmail.com>
Date: Mon, 24 Jun 2019 18:05:53 +0200
Subject: [PATCH 2/6] SVDAlgorithms -> Algorithms

---
 stdlib/LinearAlgebra/src/LinearAlgebra.jl |  6 ++++++
 stdlib/LinearAlgebra/src/svd.jl           | 25 +++++++++--------------
 stdlib/LinearAlgebra/test/svd.jl          |  2 +-
 3 files changed, 17 insertions(+), 16 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index d73807e83a0ae..a5f9136eb7b3d 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -156,6 +156,12 @@ else
     const BlasInt = Int32
 end
 
+
+abstract type Algorithm end
+struct DivideAndConquer <: Algorithm end
+struct GolubReinsch <: Algorithm end
+
+
 # Check that stride of matrix/vector is 1
 # Writing like this to avoid splatting penalty when called with multiple arguments,
 # see PR 16416
diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 7a33e31d57cbc..22034c4db03b6 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -55,11 +55,6 @@ function SVD{T}(U::AbstractArray, S::AbstractVector{Tr}, Vt::AbstractArray) wher
 end
 
 
-abstract type SVDAlgorithm end
-struct DivideAndConquer <: SVDAlgorithm end
-struct Simple <: SVDAlgorithm end
-
-
 # iteration for destructuring into components
 Base.iterate(S::SVD) = (S.U, Val(:S))
 Base.iterate(S::SVD, ::Val{:S}) = (S.S, Val(:V))
@@ -67,7 +62,7 @@ Base.iterate(S::SVD, ::Val{:V}) = (S.V, Val(:done))
 Base.iterate(S::SVD, ::Val{:done}) = nothing
 
 """
-    svd!(A; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) -> SVD
+    svd!(A; full::Bool = false, alg::Algorithm = DivideAndConquer()) -> SVD
 
 `svd!` is the same as [`svd`](@ref), but saves space by
 overwriting the input `A`, instead of creating a copy.
@@ -98,14 +93,14 @@ julia> A
   0.0       0.0  -2.0  0.0  0.0
 ```
 """
-function svd!(A::StridedMatrix{T}; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) where T<:BlasFloat
+function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = DivideAndConquer()) where T<:BlasFloat
     m,n = size(A)
     if m == 0 || n == 0
         u,s,vt = (Matrix{T}(I, m, full ? m : n), real(zeros(T,0)), Matrix{T}(I, n, n))
     else
         if typeof(alg) == DivideAndConquer
             u,s,vt = LAPACK.gesdd!(full ? 'A' : 'S', A)
-        elseif typeof(alg) == Simple
+        elseif typeof(alg) == GolubReinsch
             c = full ? 'A' : 'S'
             u,s,vt = LAPACK.gesvd!(c, c, A)
         else
@@ -116,7 +111,7 @@ function svd!(A::StridedMatrix{T}; full::Bool = false, alg::SVDAlgorithm = Divid
 end
 
 """
-    svd(A; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) -> SVD
+    svd(A; full::Bool = false, alg::Algorithm = DivideAndConquer()) -> SVD
 
 Compute the singular value decomposition (SVD) of `A` and return an `SVD` object.
 
@@ -134,7 +129,7 @@ and `V` is `N \\times N`, while in the thin factorization `U` is `M
 number of singular values.
 
 If `alg = DivideAndConquer()` (default) a divide-and-conquer algorithm is used to calculate the SVD.
-One can set `alg = Simple()` to use a simple (typically slower) algorithm instead.
+Another (typically slower) option is `alg = GolubReinsch()`.
 
 # Examples
 ```jldoctest
@@ -160,20 +155,20 @@ julia> u == F.U && s == F.S && v == F.V
 true
 ```
 """
-function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer()) where T
+function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = DivideAndConquer()) where T
     svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
 end
-function svd(x::Number; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+function svd(x::Number; full::Bool = false, alg::Algorithm = DivideAndConquer())
     SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
 end
-function svd(x::Integer; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+function svd(x::Integer; full::Bool = false, alg::Algorithm = DivideAndConquer())
     svd(float(x), full = full, alg = alg)
 end
-function svd(A::Adjoint; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+function svd(A::Adjoint; full::Bool = false, alg::Algorithm = DivideAndConquer())
     s = svd(A.parent, full = full, alg = alg)
     return SVD(s.Vt', s.S, s.U')
 end
-function svd(A::Transpose; full::Bool = false, alg::SVDAlgorithm = DivideAndConquer())
+function svd(A::Transpose; full::Bool = false, alg::Algorithm = DivideAndConquer())
     s = svd(A.parent, full = full, alg = alg)
     return SVD(transpose(s.Vt), s.S, transpose(s.U))
 end
diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl
index 0a23720bd6cc5..4ef155d7ed515 100644
--- a/stdlib/LinearAlgebra/test/svd.jl
+++ b/stdlib/LinearAlgebra/test/svd.jl
@@ -153,7 +153,7 @@ end
 
     x = [0.1 0.2; 0.3 0.4]
 
-    for alg in [LinearAlgebra.Simple(), LinearAlgebra.DivideAndConquer()]
+    for alg in [LinearAlgebra.GolubReinsch(), LinearAlgebra.DivideAndConquer()]
         sx1 = svd(x, alg = alg)
         @test sx1.U * Diagonal(sx1.S) * sx1.Vt ≊ x
         @test sx1.V * sx1.Vt ≊ I

From 8b600b7039ae290b8fc024473f67e3941dc7c990 Mon Sep 17 00:00:00 2001
From: Carsten Bauer <crstnbr@gmail.com>
Date: Mon, 12 Aug 2019 17:45:50 +0200
Subject: [PATCH 3/6] default_svd_alg

---
 stdlib/LinearAlgebra/src/svd.jl  | 22 +++++++++++++---------
 stdlib/LinearAlgebra/test/svd.jl | 11 ++---------
 2 files changed, 15 insertions(+), 18 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 22034c4db03b6..7448c9f6dff36 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -61,8 +61,12 @@ Base.iterate(S::SVD, ::Val{:S}) = (S.S, Val(:V))
 Base.iterate(S::SVD, ::Val{:V}) = (S.V, Val(:done))
 Base.iterate(S::SVD, ::Val{:done}) = nothing
 
+
+default_svd_alg(A) = DivideAndConquer()
+
+
 """
-    svd!(A; full::Bool = false, alg::Algorithm = DivideAndConquer()) -> SVD
+    svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD
 
 `svd!` is the same as [`svd`](@ref), but saves space by
 overwriting the input `A`, instead of creating a copy.
@@ -93,7 +97,7 @@ julia> A
   0.0       0.0  -2.0  0.0  0.0
 ```
 """
-function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = DivideAndConquer()) where T<:BlasFloat
+function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where T<:BlasFloat
     m,n = size(A)
     if m == 0 || n == 0
         u,s,vt = (Matrix{T}(I, m, full ? m : n), real(zeros(T,0)), Matrix{T}(I, n, n))
@@ -111,7 +115,7 @@ function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = DivideAn
 end
 
 """
-    svd(A; full::Bool = false, alg::Algorithm = DivideAndConquer()) -> SVD
+    svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD
 
 Compute the singular value decomposition (SVD) of `A` and return an `SVD` object.
 
@@ -128,7 +132,7 @@ and `V` is `N \\times N`, while in the thin factorization `U` is `M
 \\times K` and `V` is `N \\times K`, where `K = \\min(M,N)` is the
 number of singular values.
 
-If `alg = DivideAndConquer()` (default) a divide-and-conquer algorithm is used to calculate the SVD.
+If `alg = DivideAndConquer()` a divide-and-conquer algorithm is used to calculate the SVD.
 Another (typically slower) option is `alg = GolubReinsch()`.
 
 # Examples
@@ -155,20 +159,20 @@ julia> u == F.U && s == F.S && v == F.V
 true
 ```
 """
-function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = DivideAndConquer()) where T
+function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where T
     svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
 end
-function svd(x::Number; full::Bool = false, alg::Algorithm = DivideAndConquer())
+function svd(x::Number; full::Bool = false, alg::Algorithm = default_svd_alg(x))
     SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
 end
-function svd(x::Integer; full::Bool = false, alg::Algorithm = DivideAndConquer())
+function svd(x::Integer; full::Bool = false, alg::Algorithm = default_svd_alg(x))
     svd(float(x), full = full, alg = alg)
 end
-function svd(A::Adjoint; full::Bool = false, alg::Algorithm = DivideAndConquer())
+function svd(A::Adjoint; full::Bool = false, alg::Algorithm = default_svd_alg(A))
     s = svd(A.parent, full = full, alg = alg)
     return SVD(s.Vt', s.S, s.U')
 end
-function svd(A::Transpose; full::Bool = false, alg::Algorithm = DivideAndConquer())
+function svd(A::Transpose; full::Bool = false, alg::Algorithm = default_svd_alg(A))
     s = svd(A.parent, full = full, alg = alg)
     return SVD(transpose(s.Vt), s.S, transpose(s.U))
 end
diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl
index 4ef155d7ed515..c2b6b666c6b14 100644
--- a/stdlib/LinearAlgebra/test/svd.jl
+++ b/stdlib/LinearAlgebra/test/svd.jl
@@ -144,13 +144,6 @@ end
 @testset "SVD Algorithms" begin
     ≊(x,y) = isapprox(x,y,rtol=1e-15)
 
-    allpos = (v) -> begin
-        for e in v
-            e < 0 && return false
-        end
-        return true
-    end
-
     x = [0.1 0.2; 0.3 0.4]
 
     for alg in [LinearAlgebra.GolubReinsch(), LinearAlgebra.DivideAndConquer()]
@@ -158,13 +151,13 @@ end
         @test sx1.U * Diagonal(sx1.S) * sx1.Vt ≊ x
         @test sx1.V * sx1.Vt ≊ I
         @test sx1.U * sx1.U' ≊ I
-        @test allpos(sx1.S)
+        @test all(sx1.S .≥ 0)
 
         sx2 = svd!(copy(x), alg = alg)
         @test sx2.U * Diagonal(sx2.S) * sx2.Vt ≊ x
         @test sx2.V * sx2.Vt ≊ I
         @test sx2.U * sx2.U' ≊ I
-        @test allpos(sx2.S)
+        @test all(sx2.S .≥ 0)
     end
 end
 

From 44694c0ab7cf8bbe298fb8c614c1dc421787761a Mon Sep 17 00:00:00 2001
From: Carsten Bauer <crstnbr@gmail.com>
Date: Tue, 13 Aug 2019 12:29:13 +0200
Subject: [PATCH 4/6] refined docstring

Co-Authored-By: Andreas Noack <andreas@noack.dk>
---
 stdlib/LinearAlgebra/src/svd.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 7448c9f6dff36..cfca624e1a837 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -133,7 +133,7 @@ and `V` is `N \\times N`, while in the thin factorization `U` is `M
 number of singular values.
 
 If `alg = DivideAndConquer()` a divide-and-conquer algorithm is used to calculate the SVD.
-Another (typically slower) option is `alg = GolubReinsch()`.
+Another (typically slower but more accurate) option is `alg = GolubReinsch()`.
 
 # Examples
 ```jldoctest

From 6ba477887956ae51498ed8cc8faf117b73ece5c8 Mon Sep 17 00:00:00 2001
From: Carsten Bauer <crstnbr@gmail.com>
Date: Tue, 13 Aug 2019 13:18:15 +0200
Subject: [PATCH 5/6] rename to QRIteration; _svd! dispatch

---
 stdlib/LinearAlgebra/src/LinearAlgebra.jl |  2 +-
 stdlib/LinearAlgebra/src/svd.jl           | 21 ++++++++++++---------
 stdlib/LinearAlgebra/test/svd.jl          |  2 +-
 3 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index a5f9136eb7b3d..4ffb18a4d021a 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -159,7 +159,7 @@ end
 
 abstract type Algorithm end
 struct DivideAndConquer <: Algorithm end
-struct GolubReinsch <: Algorithm end
+struct QRIteration <: Algorithm end
 
 
 # Check that stride of matrix/vector is 1
diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 7448c9f6dff36..6a948807ec673 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -102,18 +102,21 @@ function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = default_
     if m == 0 || n == 0
         u,s,vt = (Matrix{T}(I, m, full ? m : n), real(zeros(T,0)), Matrix{T}(I, n, n))
     else
-        if typeof(alg) == DivideAndConquer
-            u,s,vt = LAPACK.gesdd!(full ? 'A' : 'S', A)
-        elseif typeof(alg) == GolubReinsch
-            c = full ? 'A' : 'S'
-            u,s,vt = LAPACK.gesvd!(c, c, A)
-        else
-            throw(ArgumentError("Unsupported value for `alg` keyword."))
-        end
+        u,s,vt = _svd!(A,full,alg)
     end
     SVD(u,s,vt)
 end
 
+
+_svd!(A::StridedMatrix{T}, full::Bool, alg::Algorithm) where T<:BlasFloat = throw(ArgumentError("Unsupported value for `alg` keyword."))
+_svd!(A::StridedMatrix{T}, full::Bool, alg::DivideAndConquer) where T<:BlasFloat = LAPACK.gesdd!(full ? 'A' : 'S', A)
+function _svd!(A::StridedMatrix{T}, full::Bool, alg::QRIteration) where T<:BlasFloat
+    c = full ? 'A' : 'S'
+    u,s,vt = LAPACK.gesvd!(c, c, A)
+end
+
+
+
 """
     svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD
 
@@ -133,7 +136,7 @@ and `V` is `N \\times N`, while in the thin factorization `U` is `M
 number of singular values.
 
 If `alg = DivideAndConquer()` a divide-and-conquer algorithm is used to calculate the SVD.
-Another (typically slower) option is `alg = GolubReinsch()`.
+Another (typically slower) option is `alg = QRIteration()`.
 
 # Examples
 ```jldoctest
diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl
index c2b6b666c6b14..2e882de537da1 100644
--- a/stdlib/LinearAlgebra/test/svd.jl
+++ b/stdlib/LinearAlgebra/test/svd.jl
@@ -146,7 +146,7 @@ end
 
     x = [0.1 0.2; 0.3 0.4]
 
-    for alg in [LinearAlgebra.GolubReinsch(), LinearAlgebra.DivideAndConquer()]
+    for alg in [LinearAlgebra.QRIteration(), LinearAlgebra.DivideAndConquer()]
         sx1 = svd(x, alg = alg)
         @test sx1.U * Diagonal(sx1.S) * sx1.Vt ≊ x
         @test sx1.V * sx1.Vt ≊ I

From 6918e0a1d0150011490ed59e00b9b4723045809f Mon Sep 17 00:00:00 2001
From: Carsten Bauer <crstnbr@gmail.com>
Date: Thu, 15 Aug 2019 19:41:33 +0200
Subject: [PATCH 6/6] compat annotation

---
 stdlib/LinearAlgebra/src/svd.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 3ddc64be0e996..c650eae3551b9 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -138,6 +138,9 @@ number of singular values.
 If `alg = DivideAndConquer()` a divide-and-conquer algorithm is used to calculate the SVD.
 Another (typically slower but more accurate) option is `alg = QRIteration()`.
 
+!!! compat "Julia 1.3"
+    The `alg` keyword argument requires Julia 1.3 or later.
+
 # Examples
 ```jldoctest
 julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]