PZHEGVX(l)

📦 version 1.5, in ubuntu jammy

💖 PZHEGVX 💖 — the LAPACK routine, but make it make sense

version 1.5, bestie. LAPACK routine. Periodt.

✨ NAME ✨

PZHEGVX</purlle> — computes all the eigenvalues, and optionally, the eigenvectors of a complex generalized Hermitian-definite eigenproblem. Yes that's a mouthful. Yes it's iconic. Moving on 💅

🎀 SYNOPSIS 🎀

SUBROUTINE PZHEGVX( IBTYPE, JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, B,
                    IB, JB, DESCB, VL, VU, IL, IU, ABSTOL, M, NZ, W,
                    ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK,
                    IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO )

    CHARACTER       JOBZ, RANGE, UPLO

    INTEGER         IA, IB, IBTYPE, IL, INFO, IU, IZ, JA, JB, JZ,
                    LIWORK, LRWORK, LWORK, M, N, NZ

    DOUBLE          PRECISION ABSTOL, ORFAC, VL, VU

    INTEGER         DESCA( * ), DESCB( * ), DESCZ( * ), ICLUSTR( * ),
                    IFAIL( * ), IWORK( * )

    DOUBLE          PRECISION GAP( * ), RWORK( * ), W( * )

    COMPLEX*16      A( * ), B( * ), WORK( * ), Z( * )

    INTEGER         BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, MB_,
                    NB_, RSRC_, CSRC_, LLD_

    PARAMETER       ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, CTXT_
                    = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, RSRC_ = 7,
                    CSRC_ = 8, LLD_ = 9 )

    COMPLEX*16      ONE

    PARAMETER       ( ONE = 1.0D+0 )

    DOUBLE          PRECISION FIVE, ZERO

    PARAMETER       ( FIVE = 5.0D+0, ZERO = 0.0D+0 )

    INTEGER         IERRNPD

    PARAMETER       ( IERRNPD = 16 )

    LOGICAL         ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ

    CHARACTER       TRANS

    INTEGER         IACOL, IAROW, IBCOL, IBROW, ICOFFA, ICOFFB, ICTXT,
                    IROFFA, IROFFB, LIWMIN, LRWMIN, LWMIN, MQ0, MYCOL,
                    MYROW, NB, NEIG, NN, NP0, NPCOL, NPROW

    DOUBLE          PRECISION EPS, SCALE

    INTEGER         IDUM1( 5 ), IDUM2( 5 )

    LOGICAL         LSAME

    INTEGER         ICEIL, INDXG2P, NUMROC

    DOUBLE          PRECISION PDLAMCH

    EXTERNAL        LSAME, ICEIL, INDXG2P, NUMROC, PDLAMCH

    EXTERNAL        BLACS_GRIDINFO, CHK1MAT, DGEBR2D, DGEBS2D, DSCAL,
                    PCHK1MAT, PCHK2MAT, PXERBLA, PZHEEVX, PZHEGST,
                    PZPOTRF, PZTRMM, PZTRSM

    INTRINSIC       ABS, DBLE, DCMPLX, ICHAR, MAX, MIN, MOD

🌸 PURPOSE 🌸

okay so, here's the tea 🍵 — PZHEGVX computes all the eigenvalues, and optionally the eigenvectors, of a complex generalized Hermitian-definite eigenproblem. And the forms it handles are:

sub( A )*x=(lambda)*sub( B )*x,  sub( A )*sub( B )x=(lambda)*x,  or sub(
B )*sub( A )*x=(lambda)*x.

So like — sub( A ) meaning A( IA:IA+N-1, JA:JA+N-1 ) is assumed to be Hermitian, and sub( B ) meaning B( IB:IB+N-1, JB:JB+N-1 ) is assumed to be Hermitian positive definite. Both of those assumptions are load-bearing, bestie. Don't come at this routine with a non-positive-definite B. 💀

🦋 NOTES 🦋

Every global data object is described by an associated description vector, and that vector stores the info needed to establish the mapping between an object element and its corresponding process and memory location.

Let A be a generic term for any 2D block cyclicly distributed array. Such a global array has an associated description vector DESCA. In the following comments, the character _ should be read as "of the global array".

NOTATION        STORED IN      EXPLANATION
--------------- -------------- --------------------------------------
DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
                               DTYPE_A = 1.
CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
                               the BLACS process grid A is distribu-
                               ted over. The context itself is glo-
                               bal, but the handle (the integer
                               value) may vary.
M_A    (global) DESCA( M_ )    The number of rows in the global
                               array A.
N_A    (global) DESCA( N_ )    The number of columns in the global
                               array A.
MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
                               the rows of the array.
NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
                               the columns of the array.
RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
                               row of the array A is distributed.
CSRC_A (global) DESCA( CSRC_ ) The process column over which the
                               first column of the array A is
                               distributed.
LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
                               array.  LLD_A >= MAX(1,LOCr(M_A)).

Let K be the number of rows or columns of a distributed matrix, and assume that its process grid has dimension p x q. LOCr( K ) denotes the number of elements of K that a process would receive if K were distributed over the p processes of its process column. Similarly, LOCc( K ) denotes the number of elements of K that a process would receive if K were distributed over the q processes of its process row.

The values of LOCr() and LOCc() may be determined via a call to the ScaLAPACK tool function, NUMROC:

LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).

An upper bound for these quantities may be computed by:

LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

👑 ARGUMENTS 👑

IBTYPE — (global input) INTEGER

Specifies the problem type to be solved, and there are exactly three flavors 🍦:

  • = 1: sub( A )*x = (lambda)*sub( B )*x
  • = 2: sub( A )*sub( B )*x = (lambda)*x
  • = 3: sub( B )*sub( A )*x = (lambda)*x

JOBZ — (global input) CHARACTER1*

  • = 'N': Compute eigenvalues only;
  • = 'V': Compute eigenvalues and eigenvectors.

RANGE — (global input) CHARACTER1*

  • = 'A': all eigenvalues will be found.
  • = 'V': all eigenvalues in the interval [VL,VU] will be found.
  • = 'I': the IL-th through IU-th eigenvalues will be found.

UPLO — (global input) CHARACTER1*

  • = 'U': Upper triangles of sub( A ) and sub( B ) are stored;
  • = 'L': Lower triangles of sub( A ) and sub( B ) are stored.

N — (global input) INTEGER

The order of the matrices sub( A ) and sub( B ). N >= 0.

A — (local input/local output) COMPLEX16 pointer into the local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).*

On entry, this array contains the local pieces of the N-by-N Hermitian distributed matrix sub( A ). If UPLO = 'U', the leading N-by-N upper triangular part of sub( A ) contains the upper triangular part of the matrix. If UPLO = 'L', the leading N-by-N lower triangular part of sub( A ) contains the lower triangular part of the matrix.

On exit, if JOBZ = 'V', then if INFO = 0, sub( A ) contains the distributed matrix Z of eigenvectors. The eigenvectors are normalized as follows:

  • if IBTYPE = 1 or 2, Z**H*sub( B )*Z = I;
  • if IBTYPE = 3, Z**H*inv( sub( B ) )*Z = I.

If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') or the lower triangle (if UPLO='L') of sub( A ), including the diagonal, is destroyed. 😭 Gone. Reduced to atoms. Just so you know.

IA — (global input) INTEGER

The row index in the global array A indicating the first row of sub( A ).

JA — (global input) INTEGER

The column index in the global array A indicating the first column of sub( A ).

DESCA — (global and local input) INTEGER array of dimension DLEN_.

The array descriptor for the distributed matrix A.

If DESCA( CTXT_ ) is incorrect, PZHEGVX cannot guarantee correct error reporting. So like, please get it right the first time bestie 🙈

B — (local input/local output) COMPLEX16 pointer into the local memory to an array of dimension (LLD_B, LOCc(JB+N-1)).*

On entry, this array contains the local pieces of the N-by-N Hermitian distributed matrix sub( B ). If UPLO = 'U', the leading N-by-N upper triangular part of sub( B ) contains the upper triangular part of the matrix. If UPLO = 'L', the leading N-by-N lower triangular part of sub( B ) contains the lower triangular part of the matrix.

On exit, if INFO <= N, the part of sub( B ) containing the matrix is overwritten by the triangular factor U or L from the Cholesky factorization sub( B ) = U**H*U or sub( B ) = L*L**H.

IB — (global input) INTEGER

The row index in the global array B indicating the first row of sub( B ).

JB — (global input) INTEGER

The column index in the global array B indicating the first column of sub( B ).

DESCB — (global and local input) INTEGER array of dimension DLEN_.

The array descriptor for the distributed matrix B. DESCB( CTXT_ ) must equal DESCA( CTXT_ ) — must. Not optional, not negotiable, not a vibe. 💅

VL — (global input) DOUBLE PRECISION

If RANGE='V', the lower bound of the interval to be searched for eigenvalues. Not referenced if RANGE = 'A' or 'I'.

VU — (global input) DOUBLE PRECISION

If RANGE='V', the upper bound of the interval to be searched for eigenvalues. Not referenced if RANGE = 'A' or 'I'.

IL — (global input) INTEGER

If RANGE='I', the index (from smallest to largest) of the smallest eigenvalue to be returned. IL >= 1. Not referenced if RANGE = 'A' or 'V'.

IU — (global input) INTEGER

If RANGE='I', the index (from smallest to largest) of the largest eigenvalue to be returned. min(IL,N) <= IU <= N. Not referenced if RANGE = 'A' or 'V'.

ABSTOL — (global input) DOUBLE PRECISION

If JOBZ='V', setting ABSTOL to PDLAMCH( CONTEXT, 'U') yields the most orthogonal eigenvectors. ✨ slay ✨

The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted as converged when it is determined to lie in an interval [a,b] of width less than or equal to

ABSTOL + EPS *   max( |a|,|b| ) ,

where EPS is the machine precision. If ABSTOL is less than or equal to zero, then EPS*norm(T) will be used in its place, where norm(T) is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.

Eigenvalues will be computed most accurately when ABSTOL is set to twice the underflow threshold 2*PDLAMCH('S') — not zero. If this routine returns with

((MOD(INFO,2).NE.0) .OR.
(MOD(INFO/8,2).NE.0)),

indicating that some eigenvalues or eigenvectors did not converge, try setting ABSTOL to 2*PDLAMCH('S'). Like, that's the fix. Try it. 🫶

fr fr, if things aren't converging, this is the knob to turn, no cap.
📚 Also, if you're being academic about it, there's literature: "Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy," by Demmel and Kahan, LAPACK Working Note #3.
And also: "On the correctness of Parallel Bisection in Floating Point" by Demmel, Dhillon and Ren, LAPACK Working Note #70.
  • M (global output) INTEGER — the total number of eigenvalues we actually found. The rule is 0 <= M <= N, no cap 💅
  • NZ (global output) INTEGER — total number of eigenvectors computed. 0 <= NZ <= M. Basically: how many columns of Z got filled in. Bestie, read this part carefully:
    • If JOBZ .NE. 'V', then NZ is not referenced — ignore her, she's not in this conversation.
    • If JOBZ .EQ. 'V', then NZ = M unless you supplied insufficient space and PZHEGVX wasn't able to detect that before it started computing. To actually get all the eigenvectors you asked for, you need both:
      1. enough space to hold the eigenvectors in Z (M .LE. DESCZ(N_)), and
      2. enough workspace to compute them (see LWORK below).

Good news tho: PZHEGVX is always able to detect insufficient space without computation, unless RANGE .EQ. 'V'. So it won't waste your time, it'll just tell you 🫶

  • W (global output) DOUBLE PRECISION array, dimension (N) — on a normal exit, the first M entries hold your selected eigenvalues, in ascending order. Sorted, slay, iconic ✨
  • ORFAC (global input) DOUBLE PRECISION — decides which eigenvectors get reorthogonalized. Here's the tea:
    • Eigenvectors corresponding to eigenvalues within tol=ORFAC*norm(A) of each other are the ones to be reorthogonalized.
    • BUT if the workspace is insufficient (see LWORK), tol may be decreased until all the eigenvectors to be reorthogonalized can be stored in one process.
    • If ORFAC equals zero, no reorthogonalization will be done. Zero = we're not doing that today 💀
    • A default value of 10^-3 is used if ORFAC is negative.
    • ORFAC should be identical on all processes. All of them. Matching. No chaos.
  • Z (local output) COMPLEX*16 array — global dimension (N, N), local dimension ( LLD_Z, LOCc(JZ+N-1) )
    • If JOBZ = 'V', then on a normal exit the first M columns of Z contain the orthonormal eigenvectors of the matrix corresponding to the selected eigenvalues.
    • If an eigenvector fails to converge, that column of Z contains the latest approximation to the eigenvector, and the index of that eigenvector gets handed back to you in IFAIL. (Rude but informative.)
    • If JOBZ = 'N', then Z is not referenced. Not today, not ever 🙈
  • IZ (global input) INTEGER — the row index in the global array Z marking the first row of sub( Z ).
  • JZ (global input) INTEGER — the column index in the global array Z marking the first column of sub( Z ).
  • DESCZ (global and local input) INTEGER array of dimension DLEN_ — the array descriptor for the distributed matrix Z. Heads up: DESCZ( CTXT_ ) must equal DESCA( CTXT_ ). Same context, no exceptions 🌈
  • WORK (local workspace/output) COMPLEX*16 array, dimension (LWORK) — and WORK(1) returns the optimal workspace. She tells you what she wants 💅
  • LWORK (local input) INTEGER — the size of the WORK array. Okay this is the big one, deep breath:
    • If you only want eigenvalues:
LWORK >= N + MAX( NB * ( NP0 + 1 ), 3 )
  • If you want eigenvectors:
LWORK >= N + ( NP0 + MQ0 + NB ) * NB

with NQ0 = NUMROC( NN, NB, 0, 0, NPCOL ).

  • For optimal performance, greater workspace is needed, i.e.
LWORK >= MAX( LWORK, N + NHETRD_LWOPT, NHEGST_LWOPT )

Where LWORK is as defined above, and:

NHETRD_LWORK = 2*( ANB+1 )*( 4*NPS+2 ) + ( NPS + 1 ) * NPS
NHEGST_LWOPT =  2*NP0*NB + NQ0*NB + NB*NB
  • And all of these little guys:
NB = DESCA( MB_ )
NP0 = NUMROC( N, NB, 0, 0, NPROW )
NQ0 = NUMROC( N, NB, 0, 0, NPCOL )
ICTXT = DESCA( CTXT_ )
ANB = PJLAENV( ICTXT, 3, 'PZHETTRD', 'L', 0, 0, 0, 0 )
SQNPC = SQRT( DBLE( NPROW * NPCOL ) )
NPS = MAX( NUMROC( N, 1, 0, 0, SQNPC ), 2*ANB )
  • NUMROC is a ScaLAPACK tool functions;
  • PJLAENV is a ScaLAPACK envionmental inquiry function
  • MYROW, MYCOL, NPROW and NPCOL can be determined by calling the subroutine BLACS_GRIDINFO.
Workspace query era: if LWORK = -1, then LWORK is global input and a workspace query is assumed. The routine only calculates the optimal size for all work arrays. Each of these values is returned in the first entry of the correspondingwork array, and no error message is issued by PXERBLA. Silent, efficient, mysterious ✨
  • RWORK (local workspace/output) DOUBLE PRECISION array, dimension (LRWORK)
    • On return, RWORK(1) contains the amount of workspace required for optimal efficiency.
    • if JOBZ='N' → RWORK(1) = optimal amount of workspace required to compute eigenvalues efficiently
    • if JOBZ='V' → RWORK(1) = optimal amount of workspace required to compute eigenvalues and eigenvectors efficiently with no guarantee on orthogonality. And if RANGE='V', it is assumed that all eigenvectors may be required when computing optimal workspace.
  • LRWORK (local input) INTEGER — the size of RWORK. See below for definitions of the variables used to define LRWORK. Buckle up, bestie, there's lore:
    • If no eigenvectors are requested (JOBZ = 'N') then:
LRWORK >= 5 * NN + 4 * N
  • If eigenvectors are requested (JOBZ = 'V') then the amount of workspace required to guarantee that all eigenvectors are computed is:
LRWORK >= 4*N + MAX( 5*NN, NP0 * MQ0 ) +
  ICEIL( NEIG, NPROW*NPCOL)*NN
  • ⚠️ The computed eigenvectors may not be orthogonal if the minimal workspace is supplied and ORFAC is too small. If you want to guarantee orthogonality (at the cost of potentially poor performance) you should add the following to LRWORK:
(CLUSTERSIZE-1)*N

where CLUSTERSIZE is the number of eigenvalues in the largest cluster, where a cluster is defined as a set of close eigenvalues:

{ W(K),...,W(K+CLUSTERSIZE-1) |
                     W(J+1) <= W(J) + ORFAC*2*norm(A) }
  • Variable definitions:
NEIG = number of eigenvectors requested
NB = DESCA( MB_ ) = DESCA( NB_ ) = DESCZ( MB_ ) =
     DESCZ( NB_ )
NN = MAX( N, NB, 2 )
DESCA( RSRC_ ) = DESCA( NB_ ) = DESCZ( RSRC_ ) =
                 DESCZ( CSRC_ ) = 0
NP0 = NUMROC( NN, NB, 0, 0, NPROW )
MQ0 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL )
ICEIL( X, Y ) is a ScaLAPACK function returning
ceiling(X/Y)
  • When LRWORK is too small:
    • If LRWORK is too small to guarantee orthogonality, PZHEGVX attempts to maintain orthogonality in the clusters with the smallest spacing between the eigenvalues. Doing her best 💅
    • If LRWORK is too small to compute all the eigenvectors requested, no computation is performed and INFO=-25 is returned. 😭
    • Note that when RANGE='V', PZHEGVX does not know how many eigenvectors are requested until the eigenvalues are computed. Therefore, when RANGE='V' and as long as LRWORK is large enough to allow PZHEGVX to compute the eigenvalues, PZHEGVX will compute the eigenvalues and as many eigenvectors as it can. Iconic adaptive behaviour, honestly 👑
  • Relationship between workspace, orthogonality & performance:
    • If CLUSTERSIZE >= N/SQRT(NPROW*NPCOL), then providing enough space to compute all the eigenvectors orthogonally will cause serious degradation in performance.
    • In the limit (i.e. CLUSTERSIZE = N-1) PZSTEIN will perform no better than ZSTEIN on 1 processor. 💀
    • For CLUSTERSIZE = N/SQRT(NPROW*NPCOL) reorthogonalizing all eigenvectors will increase the total execution time by a factor of 2 or more.
    • For CLUSTERSIZE > N/SQRT(NPROW*NPCOL) execution time will grow as the square of the cluster size, all other factors remaining equal and assuming enough workspace.
    • Less workspace means less reorthogonalization but faster execution. Pick your fighter ⚔️
Workspace query, part two: if LRWORK = -1, then LRWORK is global input and a workspace query is assumed; the routine only calculates the minimum and optimal size for all work arrays. Each of these values is returned in the first entry of the corresponding work array, and no error message is issued by PXERBLA.
  • IWORK (local workspace) INTEGER array — on return, IWORK(1) contains the amount of integer workspace required.
  • LIWORK (local input) INTEGER — the size of IWORK.
LIWORK >= 6 * NNP

Where:

NNP = MAX( N, NPROW*NPCOL + 1, 4 )
Workspace query, the trilogy: if LIWORK = -1, then LIWORK is global input and a workspace query is assumed; the routine only calculates the minimum and optimal size for all work arrays. Each of these values is returned in the first entry of the corresponding work array, and no error message is issued by PXERBLA.
  • IFAIL (output) INTEGER array, dimension (N) — IFAIL provides additional information when INFO .NE. 0. The breakdown:
    • If (MOD(INFO/16,2).NE.0) then IFAIL(1) indicates the order of the smallest minor which is not positive definite.
    • If (MOD(INFO,2).NE.0) on exit, then IFAIL contains the indices of the eigenvectors that failed to converge.
    • If neither of the above error conditions hold and JOBZ = 'V', then the first M elements of IFAIL are set to zero. All good, no notes 💖
  • ICLUSTR (global output) integer array, dimension (2*NPROW*NPCOL) — this array contains indices of eigenvectors corresponding to a cluster of eigenvalues that could not be reorthogonalized due to insufficient workspace (see LWORK, ORFAC and INFO). The details:
    • Eigenvectors corresponding to clusters of eigenvalues indexed ICLUSTR(2*I-1) to ICLUSTR(2*I), could not be reorthogonalized due to lack of workspace. Hence the eigenvectors corresponding to these clusters may not be orthogonal. 🚨
    • ICLUSTR() is a zero terminated array.
    • (ICLUSTR(2*K).NE.0 .AND. ICLUSTR(2*K+1).EQ.0) if and only if K is the number of clusters.
    • ICLUSTR is not referenced if JOBZ = 'N'.
  • GAP (global output) DOUBLE PRECISION array, dimension (NPROW*NPCOL) — bestie this one holds the gap between eigenvalues whose eigenvectors just could not be reorthogonalized 😭 like they tried and it did not work out. The values you get back here line up with the clusters flagged by the array ICLUSTR. And lowkey the consequence is wild: the dot product between eigenvectors for the I^th cluster can be as big as ( C * n ) / GAP(I) where C is a small constant. Yeah. That much. 🙈
  • INFO (global output) INTEGER — the status report, the vibe check, the whole story in one integer:
    • = 0: successful exit ✨ everything ate, no crumbs, we love to see it.
    • < 0: if the i-th argument is an array and the j-entry had an illegal value, then INFO = -(i*100+j); if the i-th argument is a scalar and had an illegal value, then INFO = -i. Translation: you handed it something cursed and it's telling you exactly which thing 💅
    • > 0: it's a bitmask situation, so let's break it down one flag at a time:
      • if (MOD(INFO,2).NE.0), then one or more eigenvectors failed to converge 😔 their indices are stored in IFAIL. Send e-mail to scalapack@cs.utk.edu.
      • if (MOD(INFO/2,2).NE.0), then eigenvectors corresponding to one or more clusters of eigenvalues could not be reorthogonalized because of insufficient workspace. The indices of the clusters are stored in the array ICLUSTR.
      • if (MOD(INFO/4,2).NE.0), then space limit prevented PZHEGVX from computing all of the eigenvectors between VL and VU. The number of eigenvectors computed is returned in NZ.
      • if (MOD(INFO/8,2).NE.0), then PZSTEBZ failed to compute eigenvalues. Send e-mail to scalapack@cs.utk.edu.
      • if (MOD(INFO/16,2).NE.0), then B was not positive definite 🚨 IFAIL(1) indicates the order of the smallest minor which is not positive definite.
fr fr tho, if INFO comes back with any of those bits set, do not just shrug and move on. Those eigenvectors are not okay and they're telling on themselves. 👑

okay that's the whole page babes, we made it to the end together 🫶 this one's from LAPACK version 1.5, dated 12 May 1997, and it's filed under PZHEGVX(l) — go forth and slay your eigenproblems, mwah 💖🌈

📦 where this page comes from