I think Julia handles matrices with complex elements correctly.
My task is to modify the spectrum of a Hermitian matrix H and return just the matrix with modified spectrum. i.e. I have a function f(real_vec)->real_vec
that modifies the spectrum s(H)
of a hermitian matrix H=U[s(H)]U'
. I need the result f(H) = U[f(s(H))]U'
. I think it is possible to optimize by not computing explicitly the eigfact(H)
.
Therefore I tried to write my own eigmodif
based on the Julia realization of eigfact
. That was difficult because I was lost on the line 4816 in lapack.jl, where syevr()
is wrapped up.
I need to understand where and how, Julia has converted a COMPLEX HERMITIAN matrix to a REAL-SYMMETRIC one. Theoretically it is possible, since we have a 2n by 2n matrix J that squares to minus identity; for any n by n COMPLEX HERMITIAN matrix H we then turn it into real(H).I + imag(H).J
, or in a block form
[ real(H) -imag(H) ]
[ imag(H) real(H) ]
But how does Julia do this?
Not an expert on LAPACK, but perhaps the use of macros in the definition of the eigensolvers is unclear. From linalg/lapack.jl
(around line 4900):
# Hermitian eigensolvers
for (syev, syevr, sygvd, elty, relty) in
((:zheev_,:zheevr_,:zhegvd_,:Complex128,:Float64),
(:cheev_,:cheevr_,:chegvd_,:Complex64,:Float32))
@eval begin
# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
# * .. Scalar Arguments ..
# CHARACTER JOBZ, UPLO
⋮
⋮
So the macro code uses $syevr
as a placeholder to refer to :zheevr_
and :cheevr_
in the two passes through the loop, defining the same syevr!
for different type signatures. These are LAPACK functions dedicated to Hermitian matrices and accept complex inputs. So the meat of the calculation and complex number handling goes on inside LAPACK.