Julia wrapper for LAPACK Jacobi SVD routines

svj.jl

module SVJ
#=
Wrapper for Jacobi SVD subroutines from LAPACK
*WARNING* This is a draft. Only `joba = 'G'` has been tested.

Boilerplate copied from Julia LinearAlgebra standard library
copyright (c) 2009-2019: the Julia developers https://julialang.org

License: MIT Expat

> Permission is hereby granted, free of charge, to any person obtaining a copy
> of this software and associated documentation files (the "Software"), to deal
> in the Software without restriction, including without limitation the rights
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> copies of the Software, and to permit persons to whom the Software is
> furnished to do so, subject to the following conditions:
>
> The above copyright notice and this permission notice shall be included in all
> copies or substantial portions of the Software.
>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> SOFTWARE.

=#

using LinearAlgebra
using LinearAlgebra: chkstride1, BlasInt
using LinearAlgebra.LAPACK: liblapack, chklapackerror
using LinearAlgebra.BLAS: @blasfunc
using Base: has_offset_axes

for (gesvj, elty, relty) in
    ((:dgesvj_, :Float64, :Float64),
     (:sgesvj_, :Float32, :Float32),
     (:zgesvj_, :ComplexF64, :Float64),
     (:cgesvj_, :ComplexF32, :Float32))
    @eval begin
#       SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
#                          LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
#
#       .. Scalar Arguments ..
#       INTEGER            INFO, LDA, LDV, LWORK, LRWORK, M, MV, N
#       CHARACTER*1        JOBA, JOBU, JOBV
#       ..
#       .. Array Arguments ..
#       COMPLEX*16         A( LDA, * ),  V( LDV, * ), CWORK( LWORK )
#       DOUBLE PRECISION   RWORK( LRWORK ),  SVA( N )
        function gesvj!(joba::AbstractChar, jobu::AbstractChar,
                        jobv::AbstractChar, A::AbstractMatrix{$elty};
                        V = similar(A, $elty, (0, size(A,2))),
                        tol=zero($relty))
            @assert !has_offset_axes(A)
            @assert !has_offset_axes(V)
            chkstride1(A)
            m, n   = size(A)
            if m < n
                throw(ArgumentError("matrix A must be tall or square"))
            end
            minmn  = min(m, n)
            S      = similar(A, $relty, minmn)
            if (jobv != 'N') && (size(V,2) != n)
                throw(DimensionMismatch("matrix V must have n columns"))
            end
            mv = 0
            if jobv == 'A'
                mv = size(V,1)
            elseif jobv ∈ ('V','J')
                #resize!(V,(n,n))
                V = similar(A, $elty, (n,n))
            end
            lwork  = max(6,m+n)
            work   = Vector{$elty}(undef, lwork)
            cmplx  = eltype(A) <: Complex
            if cmplx
                lrwork = max(6,n)
                rwork = Vector{$relty}(undef, lrwork)
            end
            info   = Ref{BlasInt}()
            if cmplx
                ccall((@blasfunc($gesvj), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobA/U/V
                       Ref{BlasInt}, Ref{BlasInt}, # m,n
                       Ptr{$elty}, Ref{BlasInt}, # A,lda
                       Ptr{$relty}, Ref{BlasInt}, # SVA, mv
                       Ptr{$elty}, Ref{BlasInt}, # V, ldv
                       Ptr{$elty}, Ref{BlasInt}, # cwork, lwork
                       Ptr{$relty}, Ref{BlasInt}, # rwork, lrwork
                       Ptr{BlasInt}), # info
                      joba, jobu, jobv, m, n, A, max(1,stride(A,2)), S, mv, V, max(1,stride(V,2)),
                      work, lwork, rwork, lrwork, info)
            else
                ccall((@blasfunc($gesvj), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobA/U/V
                       Ref{BlasInt}, Ref{BlasInt}, # m,n
                       Ptr{$elty}, Ref{BlasInt}, # A,lda
                       Ptr{$relty}, Ref{BlasInt}, # SVA, mv
                       Ptr{$elty}, Ref{BlasInt}, # V, ldv
                       Ptr{$elty}, Ref{BlasInt}, # work, lwork
                       Ptr{BlasInt}), # info
                      joba, jobu, jobv, m, n, A, max(1,stride(A,2)), S, mv, V, max(1,stride(V,2)),
                      work, lwork, info)
            end
            chklapackerror(info[])

            if cmplx
                scale = rwork[1]
            else
                scale = work[1]
            end
            if scale != one($relty)
                # FIXME: defeats the whole purpose
                lmul!(scale, S)
            end
            null = similar(A, $elty, (0,0))
            if jobu ∈ ('U','F','C')
                if jobv ∈ ('V','J','A')
                    return (A, S, V)
                else
                    return (A, S, null)
                end
            elseif jobv ∈ ('V','J','A')
                return (null, S, V)
            else
                return (null, S, null)
            end
        end # function
    end # eval
end # loop over types

"""
    gesvj!(joba, jobu, jobv, A; Vinit = similar(A, (0,n))) -> (U, S, V)

Finds the singular value decomposition of `A`, `A = U * S * V'`.
See the LAPACK documentation for `dgesvj` etc.
"""
function gesvj! end

end # module