3
$\begingroup$

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?

$\endgroup$
  • 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 Vincent Jun 19 '18 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$ – Kirill Jun 21 '18 at 8:59
1
$\begingroup$

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.

$\endgroup$
  • $\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 Bush Jun 21 '18 at 10:51

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.