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?
2
-
1$\begingroup$ 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. $\endgroup$– Picaud VincentJun 19, 2018 at 22:00
-
$\begingroup$ Are you asking how best to evaluate the expression, or are you asking how best to evaluate it specifically with BLAS? $\endgroup$– KirillJun 21, 2018 at 8:59
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
1
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.
-
$\begingroup$ 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. $\endgroup$– Ian BushJun 21, 2018 at 10:51