Search code examples
julialapack

Internal manipulation of complex hermitian matrix / explain the use of "RealHermSymComplexHerm" in symmetric.jl


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?


Solution

  • 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.