Implementation for matrix multiplication with `adjoint`

General Usage Performance
1

Hi,

As I understand it, adjoint creates a lazy wrapper for an adjoint of an array, and therefore does not create a copy. I’m wondering how multiplication is implemented in this case i.e. what kind of LAPACK functionalities are called to perform e.g. adjoint(T) * T for T::Array{Float64, 2}. I looked through matmul.jl and adjtrans.jl but could not quite figure out where the implementation was.

More generally, how are lazy wrappers like Transpose, PermutedDimsArray, and view of an array handled when it comes to these numerical Linear Algebra operations?

Thanks.

2

Dense matrix multiplication is handled by the gemm routine in BLAS. Here is a link to the MKL documentation of gemm: Intel | Data Center Solutions, IoT, and PC Innovation (it won’t be that different for other implementations). As you can see it accepts arguments that let you indicate to use the transpose or adjoint of either of the input matrices. The lazy transpose and adjoint wrappers just call into BLAS with the appropriate arguments.

2 Likes
3

Thanks! That makes sense. Do you know where in Julia Base I can find the function calls to BLAS?

4

The actual calls into BLAS are in stdlib/LinearAlgebra/src/blas.jl. The definitions that are relevant for making transpose and adjoint work right or stdlib/LinearAlgebra/src/matmul.jl. There are a bunch of definitions in there that mix and matrix normal, trnapose and adjoint arrays. They call gemm_wrapper and pass ‘N’, ‘T’, or ‘C’ to indicate the right way to treat the matrices.

1 Like
5

Thanks! I will take a look.