How can I get the CS decomposition of a matrix?

Deniz · February 12, 2020, 1:22am

There seems to be no function in any well-known Julia library to obtain the CS decomposition of a matrix. There is a function CUNCSD in LAPACK, but it is not available in Julia by default. Maybe I can call that function from Julia?

longemen3000 · February 12, 2020, 2:38am

lapack comes with julia in the module LinearAlgebra.LAPACK, but i can’t find a julia interface to to CUNCSD in particular. in this file are the c calls to the lapack routines:

github.com

JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/lapack.jl

# This file is a part of Julia. License is MIT: https://julialang.org/license

module LAPACK
@doc """
Interfaces to LAPACK subroutines.
""" LAPACK

using ..LinearAlgebra.BLAS: @blasfunc, chkuplo

using ..LinearAlgebra: libblastrampoline, BlasFloat, BlasInt, LAPACKException, DimensionMismatch,
    SingularException, PosDefException, chkstride1, checksquare,triu, tril, dot

using Base: iszero, require_one_based_indexing

#Generic LAPACK error handlers
"""
Handle only negative LAPACK error codes

*NOTE* use only if the positive error code is useful.
"""

This file has been truncated. show original

if you know the conventions of calling directly a lapack routine, maybe you can a direct c call.

Also, this merits an issue on the Julia repo

Deniz · February 12, 2020, 2:49am

Why is it actually that Julia doesn’t have an interface for all of LAPACK? The wrappers must be automatically generated, right?

stevengj · February 12, 2020, 2:51am

No, they are hand-written. Writing high-level wrappers for a Fortran library like this is non-trivial to automate, especially since there are four variants of every subroutine (real/complex and single/double precision), sometimes with subtly different calling conventions (e.g. different workspace sizes).

Of course, nothing prevents you from adding a wrapper for any LAPACK routine you want using ccall in user code — there’s nothing special about the code in the Julia standard library.

Deniz · February 12, 2020, 3:04am

I see. How would a wrapper for CUNCSD look like?

The specification is here: http://www.netlib.org/lapack/explore-html/d3/db9/group__complex_o_t_h_e_rcomputational_ga55ab9c7b1a2bce552037cf519201e950.html

stevengj · February 12, 2020, 3:06am

Look at the other wrappers in lapack.jl; adapting one for *uncsd shouldn’t be too hard.

dlfivefifty · February 12, 2020, 10:16am

There are a lot in LAPACK that are non-standard, for example, QL, LQ, RQ decompositions, that would be quite a bit of work to wrap.

I’ve set up MatrixFactorizations.jl as a home for other matrix factorizations. At the moment it just includes QL and a version of QR that allows more general types.