add complex bidiagonalization methods from LAPACK

This enables O(n^2) svdvals for complex matrices
JuliaLinearAlgebra · Aug 28, 2023 · acebbf1 · acebbf1
1 parent 19474cf
commit acebbf1
Show file tree

Hide file tree

Showing 4 changed files with 46 additions and 5 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "BandedMatrices"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "0.17.37"
+version = "0.17.38"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

diff --git a/src/banded/BandedMatrix.jl b/src/banded/BandedMatrix.jl
@@ -798,6 +798,20 @@ function _bidiagonalize!(A::AbstractMatrix{T}, M::BandedColumnMajor) where T
     Bidiagonal(d, e, :U)
 end
 
+function _bidiagonalize!(A::AbstractMatrix{T}, M::BandedColumnMajor) where T <: Complex
+    m, n = size(A)
+    mn = min(m, n)
+    d = Vector{real(T)}(undef, mn)
+    e = Vector{real(T)}(undef, mn-1)
+    Q = Matrix{T}(undef, 0, 0)
+    Pt = Matrix{T}(undef, 0, 0)
+    C = Matrix{T}(undef, 0, 0)
+    work = Vector{T}(undef, 2*max(m, n))
+    rwork = Vector{real(T)}(undef, 2*max(m, n))
+    gbbrd!('N', m, n, 0, bandwidth(A, 1), bandwidth(A, 2), bandeddata(A), d, e, Q, Pt, C, work, rwork)
+    Bidiagonal(d, e, :U)
+end
+
 bidiagonalize!(A::AbstractMatrix) = _bidiagonalize!(A, MemoryLayout(typeof(A)))
 bidiagonalize(A::AbstractMatrix) = bidiagonalize!(copy(A))
 

diff --git a/src/lapack.jl b/src/lapack.jl
@@ -190,6 +190,31 @@ for (fname, elty) in ((:dgbbrd_,:Float64),
     end
 end
 
+for (fname, elty) in ((:zgbbrd_,:ComplexF64),
+                      (:cgbbrd_,:ComplexF32))
+    @eval begin
+        local Relty = real($elty)
+        function gbbrd!(vect::Char, m::Int, n::Int, ncc::Int,
+                        kl::Int, ku::Int, ab::AbstractMatrix{$elty},
+                        d::AbstractVector{Relty}, e::AbstractVector{Relty}, Q::AbstractMatrix{$elty},
+                        Pt::AbstractMatrix{$elty}, C::AbstractMatrix{$elty}, work::AbstractVector{$elty},
+                        rwork::AbstractVector{Relty})
+            info  = Ref{BlasInt}()
+            ccall((@blasfunc($fname), liblapack), Nothing,
+                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
+                 Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
+                 Ptr{Relty}, Ptr{Relty}, Ptr{$elty}, Ref{BlasInt},
+                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ptr{Relty}, Ptr{BlasInt}),
+                 vect, m, n, ncc,
+                 kl, ku, ab, max(1,stride(ab,2)),
+                 d, e, Q, max(1,stride(Q,2)),
+                 Pt, max(1,stride(Pt,2)), C, max(1,stride(C,2)), work, rwork, info)
+            LAPACK.chklapackerror(info[])
+            d, e, Q, Pt, C
+        end
+    end
+end
+
 # All the eigenvalues and, optionally, eigenvectors of a real symmetric band matrix A.
 for (fname, elty) in ((:dsbev_,:Float64),
                       (:ssbev_,:Float32))

diff --git a/test/test_banded.jl b/test/test_banded.jl
@@ -451,10 +451,12 @@ Base.similar(::MyMatrix, ::Type{T}, m::Int, n::Int) where T = MyMatrix{T}(undef,
     end
 
     @testset "induced norm" begin
-        A = brand(Float64, 12, 10, 2, 3)
-        B = Matrix(A)
-	    @test opnorm(A) ≈ opnorm(B)
-	    @test cond(A) ≈ cond(B)
+        for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
+            A = brand(T, 12, 10, 2, 3)
+            B = Matrix(A)
+    	    @test opnorm(A) ≈ opnorm(B)
+    	    @test cond(A) ≈ cond(B)
+        end
     end
 
     @testset "triu/tril" begin