Add support for complex Hermitian banded matrix eigenvalue problems

[devon-research]

At the moment, it appears there is only support for eigenvalue problems for (real) symmetric banded matrices and not for complex Hermitian banded matrices.

It seems that the most direct route which mirrors the implementation for symmetric banded matrices is to implement tridiagonalization for complex Hermitian banded matrices. (In addition to sbtrd!, there should be a method hbtrd! which tridiagonalizes complex Hermitian banded matrices.)

[devon-research]

I realized after initially thinking of this issue that sbtrd! as used in this package does not actually call the LAPACK sbtrd, but instead is reimplemented in Julia. (Why was this done?) Reimplementing hbtrd sounds like a lot of work, but perhaps one could just call the LAPACK computational routine...

[MikaelSlevinsky]

Looking at #101, it's because there's a complexity benefit. For an n x n symmetric banded matrix of bandwidth m, band reduction algorithms cost O(m n²) mathematically. That is, provided that one does not apply the O(m n²) required Givens rotations to a pre-allocated matrix for dense eigenvectors. Forming the dense eigenvectors is what LAPACK does, which costs O(n³) and is much more expensive in practice than a dense eigensolve. Once one has a symmetric tridiagonal problem, then the eigendecomposition of that also costs only O(n²).

As for Hermitian problems, it's possible that the band reduction https://github.com/JuliaMatrices/BandedMatrices.jl/blob/310699438fd4f2c1089432f96e3f07b5fd47436d/src/symbanded/bandedeigen.jl#L188-L195 (and banded congruence for Hermitian-definite problems) is already implemented by multiple dispatch, though it would probably need to be checked against the Fortran reference.

[devon-research]

Forming the dense eigenvectors is what LAPACK does, which costs O(n^3) and is much more expensive in practice than a dense eigensolve. Once one has a symmetric tridiagonal problem, then the eigendecomposition of that also costs only O(n^2).

What do you mean when you say that the banded eigensolver is "much more expensive in practice than a dense eigensolve"?

Also, does this apply if one is only interested in computing the eigenvalues (without eigenvectors)?

As for Hermitian problems, it's possible that the band reduction (and banded congruence for Hermitian-definite problems) is already implemented by multiple dispatch, though it would probably need to be checked against the Fortran reference.

I was thinking this could be a possibility, but that would be very striking...

[MikaelSlevinsky]

I mean that dsbev.f is slow compared with dsyev.f (for any bandwidth). Have a look at the summary of #101 for more info. N.B. the name changes as part of the conversation (bandedeigen usurped eigen which no longer calls dsbev).

[dlfivefifty]

We should add a comment to the code explaining why we do this so we don’t forget

[devon-research]

I was just scrolling through a diff comparison of chbtrd.f and dsbtrd.f and there are two meaningful-looking chunks that are in the former and not in the latter.

Both of them have a similar form and look like this:

*           make off-diagonal elements real and copy them to E
*
            DO 100 I = 1, N - 1
               T = AB( KD, I+1 )
               ABST = ABS( T )
               AB( KD, I+1 ) = ABST
               E( I ) = ABST
               IF( ABST.NE.ZERO ) THEN
                  T = T / ABST
               ELSE
                  T = CONE
               END IF
               IF( I.LT.N-1 )
     $            AB( KD, I+2 ) = AB( KD, I+2 )*T
               IF( WANTQ ) THEN
                  CALL CSCAL( N, CONJG( T ), Q( 1, I+1 ), 1 )
               END IF

  100       CONTINUE

and replace something like this:

*           copy off-diagonal elements to E
*
            DO 100 I = 1, N - 1
               E( I ) = AB( KD, I+1 )
  100       CONTINUE

Other than that, there are some spots setting things to their real parts (I think?), some various CONJs, a lot of D*s changed to C*s, and some differing declarations.

I have never worked closely with Fortran code and the above looks rather confusing at first glance, though I am guessing I could figure it out with some Fortran syntax tutorials and a careful side-by-side read of the Julia and Fortran versions of sbtrd. If anyone thinks they can do this in 2 minutes though please say so 😅 . @dlfivefifty @MikaelSlevinsky

[KlausC]

There exists a generic implementation of tridiagonalize, which for all real and complex types, which has been developed independant of the FORTRAN sources. Maybe it is worth to integrate that into the package? https://github.com/KlausC/MatrixAlgebra.jl/blob/master/src/tridiagonalize.jl

[dlfivefifty]

A PR would be great!

[KlausC]

Will do!