# Compute bilinear form with LAPACK

I need to compute a bilinear form for a set of left and right vectors $$w_k = \sum_{i,j} V_{ik}^*A_{ij}U_{jk},$$ where $A_{ij}\in\mathbb{R}$ and $U_{jk}, V_{ik} \in \mathbb{C}$ (Assume that all the matrix sizes make sense). This can be obtained in one line of MATLAB code $$\texttt{w=diag(V'*A*U);}$$ or $$\texttt{w=sum(conj(V).*(A*U), 1);}$$ My question is, what would be the most straight forward FORTRAN implementation of this operation using proper LAPACK/BLAS, considered that $A_{ij}$ is originally given as a real array and $U_{jk}, V_{ik}$ are stored as complex arrays?

• AFAIK type mixing (real / complex) is not supported in Blas/Lapack, I fear that you will be obliged to copy A to form a complex matrix with zero imaginary part. – Picaud Vincent Jun 19 '18 at 22:00
• Are you asking how best to evaluate the expression, or are you asking how best to evaluate it specifically with BLAS? – Kirill Jun 21 '18 at 8:59

I would use ?gemm for the matrix product * (MKL reference) and v?mul for the Hadamard product .* (MKL reference). As said before, you would have to cast everything to complex to the best of my knowledge.

Supposing you work in double precision, something along the lines of:

Integer, Parameter :: dp = kind(1.0d0)
Real(dp), Dimension(:,:), Allocatable :: V_real
Complex(dp), Dimension(:,:), Allocatable :: A, U, V
Complex(dp), Dimension(:,:), Allocatable :: C, D ! some auxilliary arrays
Complex(dp), Dimension(:), Allocatable :: w
!allocate and define everything
!convert V from real to complex:
V = Cmplx(V_real, 0._dp, kind=dp)
!start by computing C = A*U with ?gemm
Call zgemm(..., A, ..., U, ... C, ...)
!Then compute w = sum(conj(V).*C) with vzmul
V = Conjg (V)
Call vzmul (n, V, C, D)
w = sum(D,1)


You may not need the auxilliary arrays C and D if you can overwrite some of the original arrays A, U, V, and if the shape are compatible (but I assumed that you could not do that).

Of course for both routine/function there are quite some book-keeping arguments that you need to figure out using the reference (essentially to indicate the logical and physical shape of your arrays). I mentionned intel's BLAS because I am familiar with it, it is free, optimized for intel machines, and the documentation is easy to browse. Further, it features v?mul, which I am not sure exists in the (non intel) regular BLAS.

• To do the element by element multiply just use a multiply sign which in standard Fortran does precisely this - vmul is not standard BLAS, and so could cause portability problems. – Ian Bush Jun 21 '18 at 10:51