CODEWITHSUNDEEP
matrixvectormultiplicationlapack
zeus300
zeus300

Reputation: 1075

Matrix-vector-multiplication, Lost in Lapack

I would like to compute

x^T A x* (x transposed times matrix A times x conjugate complex)

using Lapack. I did not find a function for this. Or is there one? What is the fastest way of doing this?

Upvotes: 0

Views: 386

Answers (1)

jack
jack

Reputation: 1790

You could compute the operation in two separate steps:

  1. compute y <- A x*
  2. compute res <- x^T y

I implemented this in the following program

program main

  use blas95

  implicit none

  integer, parameter :: n = 2
  complex :: x(n), y(n), A(n,n), res

  complex :: cdotu

  ! init vector x, matrix A by some data
  x = [(1, 1), (-1, 2)]
  A = reshape([1, 2, 3, 4], [n, n])

  ! compute: y <- A * conjg(x)
  call cgemv('N', n, n, (1, 0), A, n, conjg(x), 1, (0, 0), y, 1)    ! standard BLAS
  call gemv(A, conjg(x), y)                                         ! BLAS95 from intel MKL

  ! compute: res <- x^T * y
  res = cdotu(n, x, 1, y, 1)                                        ! standard BLAS
  res = dotu(x, y)                                                  ! BLAS95 from intel MKL
end program

Note: I included the standard BLAS routines as well as the BLAS95 routines available through intel's MKL.

Link to intel's documentation on

Upvotes: 1

Related Questions